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Abstract 

The local field correction to the spontanous dacay rate of an impurity source atom imbedded in 
a disordered dielectric is calculated to second order in the dielectric density. The result is found 
to differ from predictions associated with both "virtual" and "real" cavity models of this decay 
process. However, if the contributions from two dielectric atoms at the same position are included, 
the virtual cavity result is reproduced. 



I. INTRODUCTION 



The problem of spontaneous emission from an atom imbedded inside a dielectric has 
attracted considerable interest 111. Most theoretical treatments of this problem follow 

n 

a macroscopic approach |2J. Based on different models of the local environment of the 
imbedded atom, they give different types of local field corrections to the spontaneous de- 
cay rate T of the impurity atom. The so-called virtual cavity model gives a decay rate, 
^virtual = V^i^Y^o, assuming that a virtual cavity surrounds the emitter, while the 
"real" cavity model gives a rate, r real = y/e (^fi) 2 ^ ' assumm g that an empty spherical 
cavity surrounds the emitter. The quantity e is the permittivity of the dielectric, which is 
connected to the microscopic polarizability a by the Lorentz-Lorenz relation e = 1 + A " 



where N is the dielectric density. Expansions for the decay rates in powers of Na yield 
r rea i = [1 + iNa + f 2 {Naf + O(Na)' 3 }T and T mrtual = [1 + \Na + g(iVa) 2 + O(Na) 3 ]T . 
To first order in Na, the real and virtual cavity models give identical results, but they differ 
in higher order. To determine the validity range of these macroscopic models, calculations 
using a somewhat more fundamental approach are needed. Several attempts at such mi- 
croscopic models invo 



approach for crystals 
disordered dielectrics 



ve (i) a polariton approach for crystals jjj, (ii) a Green's function 
and disordered dielectrics ^, and (iii) an amplitude approach for 
6|0. In the polariton method, the interaction between the vacuum 
radiation field and the crystal atoms is solved exactly; the eigenmodes of this system are the 
polaritons. The source atom then decays by radiating polaritons. This polariton calculation 
agrees with the virtual cavity result |3[. In the Green's function approach, the modification 
of the decay rate results from scattering of radiation emitted from the source atom by the 
dielectric, calculated to all orders in the dielectric density. This calculation reproduces the 
virtual cavity result with the source atom at an interstitial position and the real cavity result 
with source atom at a substitutional position in the crystal^. For disordered dielectrics, 
the Green's function method gives the virtual cavity result a|. The amplitude method rep- 
resents a direct calculation of the modification of the decay rate as a perturbation series in 

no 

Na [6J 7]. To first order in Na, the radiation emitted by the source atom is scattered back 
to the source atom by a single dielectric atom; the resultant decay rate agrees with the both 
virtual and real cavity models to first order in Na J(|. 

In this paper, the amplitude method is extended to second order by including scattering 
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events in which the radiation emitted by the source atom is scattered back to the source 
atom by a combined scattering from two dielectric atoms. It will be seen that the result 
differs from those of both the real and virtual cavity models; however, when contributions 
to the decay rate originating from scattering by two dielectric atoms located at the same 
physical point are included, the calculation reverts to the virtual cavity model. 



II. CALCULATION OF SECOND ORDER CONTRIBUTION 

The source atom located at R = 0, has a J = ground state and a J = 1 excited state, 
the frequency separation of the ground and excited state denoted by Uq. The uniformly 
distributed dielectric atoms have J = ground states and J = 1 excited states, the frequency 
separation of the ground and excited state denoted by u. At t — 0, the source atom is excited 
to the m = excited state sublevel, the dielectric atoms are all in their ground states, and 
there are no photons in the field. The process we consider is radiation emitted by the source 
atom that is scattered by two dielectric atoms back to the source atom. It is assumed that 
\uj — ujq\ /r 1 but that \uj — u \ / (uj + uj ) <C 1 [rotating wave approximation (RWA)]. 



We use a multipolar Hamiltonian 



The free part is 



j m=-l 

where a z = (|2) (2| — |1) (1|), \2) and |1) are the m = excited and J = ground state 
eigenkets of the source atom, respectively, a^\m) = Hm)^ (m\ — \g)^ (g\j is the popu- 
lation difference operator between excited state \J = l,m) and ground state \ J — 0, g) of 
dielectric atom j, and is the annihilation operator for a photon having momentum k 
and polarization A. A summation convention is used, in which any repeated symbol on the 
right hand side of an equation is summed over, unless it also appears on the left-hand side 
of the equations. 

The interaction part of the Hamiltonian is, 

eo eo 

where d and dj are the dipole operators of the source atom located at the origin and a 
dielectric atom located at position Rj respectively. The operator D is the displacement field 
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having positive frequency component 



D + (R)=z eo ^j|^ e W akAeJk . R (3) 

k,A V 

where V is the quantization volume and e k is a unit polarization vector, with 

e k = cos cos k x + cos #k sin (/> k y — sin # k z (4) 
e k = - sin k x + cos k y. (5) 

In the RWA, one can write 

y = ^% k (a + a k -aj c a_)+ £ ^( 5kA (m)af (m)a kA e*- R + 5kA (m)*4 A ^' ) (m) e - ik - R ) (6) 

k k,A,m 



where the cr-t are raising and lowering operators for the source atom and a± (m) are raising 
and lowering operators between the excited state \ J — 1, m) and the ground state | J — 0,g) 
of dielectric atom j, // is the reduced matrix element of the dipole operator do and // is 

(A)^ . _ _ ^(£%±i(£\ 



that of dj between ground and excited state manifolds. The (e k )±i = =F k ^ 
(e k A ^)o = (e^ )* are spherical components of the polarization vectors. The source atom 
interacts only with the z component of the radiation field. 

The calculation proceeds as in Ref. [oj, with the addition of terms that couple dielectric 
atoms to dielectric atoms via the radiation field. After eliminating intermediate states 
involving the radiation field, one arrives at 

k = - 7 & 2 - 7 (-)e' A *G'o, mi (R i) a;)6 m .(t) (9a) 
ul 

b mj = -l'b m] - 7 (— )e~ lAt G mfi (Rj, u> )b 2 (t) - iG mjjm > s (R,- - R s , uj)b m > a {t), (9b) 
fi 

where 7 = 2/! 2 co>q/3/ic 3 and 7' = 2/i /2 c<jQ/3frc 3 are (half) the excited state decay rate of the 
source and dielectric atoms, respectively, b 2 is the state amplitude for the source atom to 
be in state |2) = \J = l,m — 0) and all dielectric atoms in their ground states, and b m . t is 
the state amplitude for dielectric atom j to be in excited state | J — 1, m) all other atoms in 



their ground states. We have set b 2 (t — r) « b 2 (t) and b nij {t — T) ps b mj (t) on the assumption 
that jRq/c, j'Rq/c <C 1, where i?o is the sample size. The quantity G mj ^(Rj — R s ,u;) is 
a propagator for scattering from a dielectric atom in sublevel rrij at position R., to one in 
sublevel m s at position R s given by 

G mj , m ,{n,u) = A 1 /"* dT r ^e-^-)- / ^^(^U^'*, (10) 

87T 7TUr J 

while G m) o(Rj, a;) is a propagator for scattering from the source atom to a dielectric atom 
in sublevel rrij at position Rj. In what follows we ignore the difference between u Q and lu, 
consistent with the RWA. 

In order to solve Eqs. (|9al9b|) . we assume that b 2 varies slowly on the time scale 1/A. If 
b mj = 2/m, e_lA *) Eqs. (j9a|9bj) are transformed to 



h = -jb 2 - t(— )G 0)m (R j ,u} )y mj (t) 
/i 

(7' - iA)y m] = -^(—)G mfi (Rj,uj )b 2 (t) - 7 / G mi)m ^(R j - R 3 ,u )y m > a (t) 
A 4 

The formal solution for bo is 



(11a) 
(lib) 



&2 = ~lh +7( — )G ,m 3 (Rj,Wo) 

A 4 



1 



Y - iA + 7'G 



7( — ) Gr m' o(R S ,^o)&2 



;i2) 



where G is an 3iV x 3N matrix having matrix elements G m . >m > g (Rj — R s ,u;o). This can be 



expanded as a power series in Na with a 



hA 



. To second order, one finds 



-762,0 



1 + iNa e^N G 0,m j (Aj 1 u )G mjt0 (Rj, UJ ) 
— (Na) 2 (-f^) 2 Go,m j (Rj,uo)G mj ,m' 3 (~R-j — R s , Uo)G m > s! o(R s , Vo) 



The term linear in the density gives the first order local field correction 
which agrees with both the virtual and real cavity models to this order. 
We now calculate the second order correction, 

n-n 



5 7 « _ 7 



(13) 



Na |6|, 



■(Na) 



"0 \2 



Go^mj (Rj; uo)G m . m i g (Rj — R s , wo)G m ' sl o (R s ,u; ). (14) 



7 'QkN' 

The sum over Rj and R s can be converted to integrals using N J dH. In this manner, 

one finds 

r5 7 ( 2 ) ^3 
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-(Nafi^f J y"dR 2 dR 1 G 0)m ,(R 2 ,a;o)G m . imi (R,a;o)G mii o(Ri,a;o) (15) 



where R =i? 2 — The next step is to evaluate the G mhm > s (R,uj). The details of the calcu- 
lation are given in the Appendix and one obtains 



Gii 
Goo 

G-1,1 

G\fi 
G-\n 



Aizh Q {k R)Y , {R) - ^ ^-h 2 {k R)Y 2> o(R); 

47T 




Anh (k R)Y 0fi (R) + \/ —h 2 (k R)Y 2>0 (R); 



3 /8tt 

2 V T5 

3 /8tt 

2 V 15 

3 /4tt 
'2V 15 



h 2 {k R)Y 2 ^ 2 {R); 

h 2 (k R)Y 2t2 (R); 

h 2 (k R)Y 2 ^ 1 (R); 



\J^h 2 {k Q R)Y 2jl {R) 



(16a) 
(16b) 
(16c) 
(16d) 
(16e) 
(16f) 



where Yg jTO (R) is a spherical harmonic and ko = ujq/c. The remaining G m . jm rs are obtained 
using G-x,-i — Gu, Go,-i — — Gift, and G ,i = — G-i,o . The spherical Hankel functions 
of the first kind, ho(koR) and h 2 (koR), conform to the appropriate boundary conditions in 
which only outgoing scattered waves are considered. 



The calculation for 



<5 7 (2) 



is tedious, since it involves contributions from nine terms. We 



will show how to calculate one specific contribution, rnj = 1, m' s = 1, and then give the final 
results for the other components. Substituting Eqs. ()16api6e|) in Eq. (fTHjl . we find 



<M 2) (m: 
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dR 1 dR 2 / i2 (/c i2 2 )F 24 (R 2 )/i (A;o^2i)^Q ) o(R2i)^(A;oi?i)l2,-i(Ri 




dR 1 dR 2 /i 2 (A:oi2 2 )r 2 ,i(R 2 )/i 2 (A:oi? 2 i)r 2 ,o(R2i)/i2(A:oit:i)y2 > -i(Ri) 

(17) 



W 5^ 

To evaluate this, we expand the Hankel functions as j^] 

hi(k Q R 21 )Y ltm (R 21 ) = i h+h ~\-lt +m ^^{2l + l)(2/i + l)(2/ 2 + 1) 



x 



- + (-1)' e(R x -R 2 








x fc h ( /c i?> ) j h ( /c #< ) Y tl mi ( R> ) Y is 



,m2 V A,J < 



R< 



(18) 



where ("') is a 3-j symbol, ji(x) is a spherical Bessel function, and R > (R<) is the larger 
(smaller) of Ri and R 2 . When this expansion is used in Eq. (fT7|) . the angular integration 



selects out only I = 2, and m = 1, — 1, terms, such that 

^ (2)(l,l) = r d P2P i r dp^ihMhzipjh^jM (19) 

7 i Jo Jo 

with pi = koRi, p 2 = koR 2 . To evaluate the above integral, we add a convergence factor 
e~ ep2 , and eventually take the limit e — > 0. The imaginary part of the integral diverges as 
P2 - ► 0, but the real part is finite and gives the local field correction to the decay rate. The 
result is 

^>(M) = jj (jVa) , 
7 n □ 

and the corresponding results for the other terms are |10|] 
^7 (2) (-l,-l) <h (2) (M) 15 (M , 2 

W) = 26 
7 oo 

57 (2) (0,1) ^7 (2) (1,0) 5 7 ( 2 )(-l,0) 57^(0,-1) 3 , 2 



77 7 7 28 



(NaY 



5 7 (2) (-l,l) 57 (2) (1,"1) 3 /Ar , a 



7 7 56 

The total second order correction to the decay rate is 

*7 (2) = + ^7 (2) (0,0) + /7 (2) (0,1) + 2 ftM) = 71 a> 

77 7 7 7 72 

This result differs from both the virtual [||(iVa) 2 ] and real [y|(iVa) 2 ] cavity models . 

Our result can be compared with Fleischhauer's j^]. The Fourier transform of 
G mim '(R, c<Jo) is the tensor field propagator F^(q, a;) in his paper, differing only by pref- 
actors. The integral (|THj) can then be done in either coordinate or momentum space. 
The momentum space integration gives a different result than our coordinate space cal- 
culation above. This surprising discrepancy can be explained by the way we expand 
hi(koR2i)Yi tm (R,2i)- The expansion we used is valid for R 1 > R 2 or R 2 > R\, but is not 
defined for Ri = R 2 . For a well-behaved integral this will not make any difference since 
Ri = R2 contributes a set of measure zero. In the present case, however, where the dipole- 
dipole interaction between dielectric atoms diverges when one atom is on the top the other, 
i.e. when Ri = R2, the contribution from R\ = R 2 can be finite. 

It is not easy to estimate this contribution in the original form of the integral 1)1 7]). 
Instead, it proves useful to Fourier transform just one of the G m>m i in the integrand. As an 



example, we consider the integral in the second term of Eq. (|17j) 

1= f /dR 2 dR 1 / i2 (fci? 2 )F 2 , 1 (R 2 )/ i2 (A;i? 21 )F 2 ,o(R 2 i)/i 2 (A; J R 1 )y 2i _ 1 (R 1 ) (21) 



We Fourier transform h2(kR 2 i)Y2fi(R2i)e~ eR21 , using a convergence factor e~ eRil that is 
physically connected with the boundary condition of outgoing spherical waves. Carrying 
out the Fourier transform in Eq. (|21j). we find 

I = ~ f [ dR 2 dR 1 h 2 (kR 2 )Y 2 ,i(R 2 )h 2 (kR l )Y^ 1 (R 1 ) [ V \ ■ , Y 2fi (f>)e^-^ 

k 6 J J J (In) 6 k z — p z + xe 

(22) 

The angular integrations can be done by expanding e* pR2 , e _ipRl in terms of spherical 
harmonics and Bessel functions. In this manner one obtains 

1 = Ti^l^ I J n^(kR,)h 2 (k Rl ) j^ ?w JL_ (23) 

We are interested only in the contribution in the region where R\ = R 2 . This contribution 
can be isolated by integrating R 2 from Ri — a to R± + a, and then integrating the resultant 
expression over p using the method of residues. In the limit that both a and e tend to zero, 
one obtains the contribution 51 from the region R\ — R 2 as 



_i , 9? f —Witt 

Re[8I] = Re[— v^-p J dR 2 R 2 2 h 2 (kR 2 )h 2 (kR 2 )} = (24) 

(the imaginary part of 51 diverges). The contribution from the sphere R\ = R 2 is identical 
to that from R x = R 2 since all other points with R x ^ R 2 on the sphere are regular and 
contribute zero to the integral. The same calculation can be done for the first integral in 
the Eq. (|17)1. For this term, there is no contribution from the region R\ = R 2 (no delta 
function like term is found) since ha(kR 2 i) has a lower order divergence at R 2 \ = than 
does h 2 (kR 2 i). 

Including contributions of the type (J2H), we find 

^' 2) < -°> = V). 

7 3 
5 7 (2) (M) 5^(1,-1) 7 2 



7 7 48 

r_i n\ W 2 )rn _-n i 

-(NaY 



5 7 (2) (0,1) 57 (2) (1>0) 5 7 (2) (-l,0) 57^(0,-1) 1 fAT , 2 



77 7 7 12 

= ^(i.-D = _i (Wo) . (25) 

When these are summed the total = ^(A^a) 2 agrees with the virtual cavity result. 



III. DISCUSSION 



The second order contribution to the modified spontaneous emission rate of an impurity 
atom in a disordered dielectric has been calculated using a microscopic theory. Depending 
on the manner in which overlapping atoms are treated, one arrives at different results. 
If the delta function contributions at Ri = R2 are included, the virtual cavity model is 
recovered, but if such terms are excluded, neither the real nor virtual cavity model results 
are found. It seems to us somewhat of an open question at this point as to whether or 
not such contributions can be uniquely calculated once Eq. (|T2*j) is expanded in a power 
series in the density. The reason for this is that the expansion parameter is not small as 
interatomic distances tend to zero. That the expansion can lead to divergences is already 
evident if the integrations are carried out using a different set of variables [10J . From physical 
considerations, however, the decay rate does not diverge, even for interparticle spacings much 

1 /3 

less than a wavelength. Actually, dielectric atoms within a sphere of radius A (7'/ A) 1 
reradiate collectively; outside this radius, there is destructive interference resulting in some 
additional finite contribution to the decay rate. In dealing with a homogeneous dielectric, 
we have performed the ensemble average by integrating over all space assuming a constant 
density. This averaging process includes configurations where interparticle spacings are 
sufficiently small to invalidate the expansion ([13)1 . Nevertheless, the procedure has yielded 
finite results for the change in the decay parameter. 

Different experiments support both the real and virtual cavity results The source 
atom in these experiments is usually an impurity ion in a protective molecular cage. No 
experiments of this nature have been carried out with impurity atomic radiators in a dielec- 
tric that consists of a dense atomic vapor. It may be possible to use an alkali metal atom 
as the source atom and rare gas atoms as the dielectric atoms. With such a system, one 
could not make the rotating wave approximation used in this paper, but the physics is not 
changed in any substantive manner. The key feature of the alkali metal - rare gas system 
is the extremely small quenching cross sections for rare gas collisions to inelastically change 
the electronic state of the alkali atom [12]. Any quenching cross sections would appear as a 
modification of the decay rate that would mask the sought after effect. For rare gas pressures 
on the order of 100 atmospheres, we estimate that a change in the decay rate of order of 3% 
could be observed. To increase the effect it is necessary to find radiator atoms whose first 
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excited state is radiatively coupled to the ground state and dielectric atoms whose lowest 
excited state is about 0.2eV above the energy of the excited state of the radiator. In this 
limit, quenching will be negligible, but the detuning A is decreased from the alkali-rare gas 
system by a factor of 50. At the same time, it is necessary to achieve a high pressure for 
the dielectric atoms. A possible system would be Li radiators with a high density sodium 
dielectric; the energy mismatch of Li and Na is about 0.25eV, giving a correction factor to 
the lithium decay rate of 1.3 x 10~ 21 iV, where N is the sodium dielectric density in units of 
atoms/cm 3 . 
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V. APPENDIX 

In this appendix, we calculate explicitly G li i(R,o;) given in Eq. (JTUJ). The other 
G mtm i(R,u) are calculated in a similar fashion. To carry out the angular integrations, one 
expands e ikR as 

i 

e ik-R = Ar[ J2 SY* n (k)Y lm (R)ji(kR), (26) 

m=— I 

uses the fact that (e^ A) )^(e^ A) )i = |(1 + cos 2 9) = ^[2Y 00 (k) + -j E Y 20 (k)}, and the orthogo- 
nality of the spherical harmonics, to obtain 

i rt poo i 

G 1a (R,uj) = -— / dr / du k ujle~^-^ T [2Y 00 (R)UkR) - —Y 20 (R)j 2 (kR)} (27) 
V 7 ^ Jo J-oc ' V5 

The spherical Bessel function can written in terms of spherical Hankel functions as ji(kR) = 
\{hi(kR) + h* l (kR)}, transforming Eq. into 

G hl (R,u) = / dr ^e- i ^-^ T {2r o(R)[^(A:i?)+/ i o(^)]-^=Y 2 o(R)[/i 2 (A ; i?)+^(fci?)]}. 

Jo J-oc V5 

In the calculation we always make the Wigner-Weisskopf approximation. Differences be- 
tween u, ujq and Uk are neglected except they appear as exponential factors. In integrating 
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over ujk, the h^(kR) terms give a contribution proportional to S(R/c + r) while the hi(kR) 
terms give a contribution proportional to 5(R/c — r). We retain only the 5(R/c — r) con- 
tributions since they correspond to the retarded solution (outgoing spherical waves). As a 
consequence, we find 
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